(* ::Package:: *)

(************************************************************************)
(* This file was generated automatically by the Mathematica front end.  *)
(* It contains Initialization cells from a Notebook file, which         *)
(* typically will have the same name as this file except ending in      *)
(* ".nb" instead of ".m".                                               *)
(*                                                                      *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs.  Doing so is equivalent   *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end.                                                                 *)
(*                                                                      *)
(* DO NOT EDIT THIS FILE.  This entire file is regenerated              *)
(* automatically each time the parent Notebook file is saved in the     *)
(* Mathematica front end.  Any changes you make to this file will be    *)
(* overwritten.                                                         *)
(************************************************************************)



BeginPackage["W`SVM`"];


SVMTrain::usage="SVMTrain[X, Y, options] use to train a SVM,
X must be a matrix(sample must be 2d at least),Y is the label.";
SVMClassify::usage="use this function to classify";

SVMSave::usage="SVMSave[s,fn_:Null] save a svm to file";
SVMLoad::usage="SVMLoad[fn_:Null] load a saved svm from file";

SVMPlot::usage="";


(**)
KernelMatrix::usage="KernelMatrix[K,X] computes the (full) matrix with elements K[x\[LeftDoubleBracket]i\[RightDoubleBracket],x\[LeftDoubleBracket]j\[RightDoubleBracket]].";

IdentityKernel::usage="IdentityKernel[x,y] is the function x.y.";
RBFKernel::usage="RBFKernel[x,y,\[Gamma]] is the function Exp[-\[Gamma](x-y).(x-y)].";
PolynomialKernel::usage="PolynomialKernel[x,y,\[Gamma],d,c], where d is an integer, is the function (\[Gamma] x.y+c\!\(\*SuperscriptBox[\()\), \(d\)]\).";
(**)
QPSolve::usage="QPSolve[X_,Y_,K_,Cs_,\[Tau]_,method_] is a wrapper func for QP problem, refer to qpType for method options";
(*QPsolves::usage="classic qpsolve, very slow";*)

(*oSVM::usage="test";*)


(*options: type of Kernel and SVM*)
svmType::usage=" =1, separable svm; =2, C-svm; =3, \[Nu]-SVM";
kernelType::usage="=1, linear; =2, PolynomialKernel; =3, RBF; ";
qpType::usage="=1, old slow for C-svm; =2, Dlib SMO for C-svm; =3, Dlib for nu-Svm";

cost::usage="the parameter C(Cp and Cn for 2-class) of C-SVC,def: {1,1}";
nu::usage="the parameter nu of nu-SVC";
t::usage="solution tolerance parameter, def: 0.001";
(*kernel function opt*)
degree::usage="def: 3";
\[Gamma]::usage="gamma in kernel, def: 1/num_feature";
coef0::usage="def:0";


(*********************************************************************)
(*we need a structure to save a trained SVM like libSVM*)(* we can use a array-like func group to describe it,
the classify func is y=sign(Underscript[\[Sum], i]ai yi K(xi, x) + b),
symbol must be local,so string is used here.
"kernel"=K, "offset"=b, "alpha"=ai
"svi"=ai * yi, "svX"= X[[sv]], "svY"=Y[[sv]](sv==ind(ai!=0))
"opts"=opt,  "label" = different labels of classes

*)
(*how to deal with errors.:
using $Failed if possible to indicate errors.
*)


(*TODO:
how to check the options,

*)
Begin["`PP`"];


(*like libSVM, scale data to [-1,1] or sth linearly*)
SVMScale[]:=Module[{},
Return["not implemented"];];

Options[SVMTrain]={svmType->1,kernelType->1,qpType->2,t->0.0001,cost->{1,1},nu->0.5,\[Gamma]->0.1,degree->3,coef0->0};
(*wrapper SVM train func*)
SVMTrain[X_?MatrixQ,Y_List,opt___]:=Module[
{m,n,K,s},
(*no tensor here*)
{m,n}=Dimensions[X];
If[Length[Y]!=m,Print["x and y dim error"];Return[]];

SetOptions[SVMTrain,\[Gamma]->1/n];
(*which kernel*)
K=Switch[(kernelType/.{opt}/.Options[SVMTrain]),
1,IdentityKernel,
2,PolynomialKernel[#1,#2,\[Gamma],degree,coef0]&/.{opt}/.Options[SVMTrain],
3,RBFKernel[#1,#2,\[Gamma]]&/.{opt}/.Options[SVMTrain],
_,Print["kernelType error, exit."];Return[$Failed];
];

(*Print["SVM train begin"];*)
(*svm is built in the specific implementation*)
s=Switch[(svmType/.{opt}/.Options[SVMTrain]),
1,oSVM[X,Y,K,opt],
2,cSVM[X,Y,K,opt],
3,SetOptions[SVMTrain,qpType->3];nuSVM[X,Y,K,opt],
_,Print["svmType error, exit."];Return[$Failed];
];
If[s==$Failed,Print["SVM train error"];Return[$Failed]];
(*TODO : although only 2-class implimented, label and groups are found from Y*)
s["label"]=Union[Y];

(*SVM["kernel"]=K;*)
s["opts"]=opt;

s
];

(*s saved SVM; data to classify, must be a matrix*)
(*save a func in SVM is not good, so the opt counts*)
SVMClassify[s_,data_?MatrixQ]:=Module[
{K,f,r,rl,opt},

opt=s["opts"];
K=Switch[(kernelType/.{opt}/.Options[SVMTrain]),
1,IdentityKernel,
2,PolynomialKernel[#1,#2,\[Gamma],degree,coef0]&/.{opt}/.Options[SVMTrain],
3,RBFKernel[#1,#2,\[Gamma]]&/.{opt}/.Options[SVMTrain],
_,Print["kernelType error, exit."];Return[$Failed];
];
(*K=s["kernel"];*)
(*\[Alpha]=s["alpha"];*)

f[x_]:=s["svi"].Map[K[#,x]&,s["svX"]]+s["offset"];
r=Map[f,data];
(*TODO: change. NOW 2-class only,
make it VF-compatible*)
rl=Sign[r];
{r,rl}
];

(*save a SVM in a file*)
SVMSave[s_,fn_:""]:=Module[
{f},
If[fn=="",
f="..\\tmp\\"<>"svm.mx",f=fn];
(*At the moment, merely save it*)
Export[f,s]
];
SVMLoad[fn_:""]:=Module[{f},
If[fn=="",
f="..\\tmp\\"<>"svm.mx",f=fn];
Import[f]
];


(*Visualization*)
(*X is the data points ploted, Y is the label of X*)
SVMPlot[s_,X_,Y_:{}]:=Module[
{K,opt,f,dim,sv,x1,x2,x3},
opt=s["opts"];
K=Switch[(kernelType/.{opt}/.Options[SVMTrain]),
1,IdentityKernel,
2,PolynomialKernel[#1,#2,\[Gamma],degree,coef0]&/.{opt}/.Options[SVMTrain],
3,RBFKernel[#1,#2,\[Gamma]]&/.{opt}/.Options[SVMTrain],
_,Print["kernelType error, exit."];Return[$Failed]
];
(*first check dim, dim must be 2 or 3*)
dim=Length[s["svX"][[1]]];
If[dim>3||dim<2,Return[$Failed]];


If[Y!={},
   1==1;
];
(*Print[x1range,x2range,x3range];*)
   f=s["svi"].Map[K[#,{x1,x2,x3}]&,s["svX"]]+s["offset"];

    Show[
      ListPointPlot3D[{s["svX"],X},PlotStyle->{Hue[0.6],PointSize[Large]}],
    (*  ListPlot[Extract[X,Position[y,1]]],
      ListPlot[Extract[X,Position[y,-1]],PlotStyle->{GrayLevel[0.6]}],*)
      ContourPlot3D[f==0,{x1,Min[X[[All,1]]],Max[X[[All,1]]]},{x2,Min[X[[All,2]]],Max[X[[All,2]]]},{x3,Min[X[[All,3]]],Max[X[[All,3]]]}],
      ContourPlot3D[f==1,{x1,Min[X[[All,1]]],Max[X[[All,1]]]},{x2,Min[X[[All,2]]],Max[X[[All,2]]]},{x3,Min[X[[All,3]]],Max[X[[All,3]]]}(*PlotStyle->Dashing[{.01,.01}]*)],
      ContourPlot3D[f==-1,{x1,Min[X[[All,1]]],Max[X[[All,1]]]},{x2,Min[X[[All,2]]],Max[X[[All,2]]]},{x3,Min[X[[All,3]]],Max[X[[All,3]]]}(*PlotStyle->Dashing[{.01,.01}]]*)],BoxRatios->Automatic
]
];



(*There are different kernels and dfferent SVMs*)
(*kernels, dots apply to list only!!all default values are same as libSVM*)
IdentityKernel[x_List,y_List]=x.y;
(*\[Gamma]>0*)
RBFKernel[x_,y_,\[Gamma]_]=E^(-\[Gamma](x-y).(x-y));
(*d \[Element] \[DoubleStruckCapitalZ]*)
PolynomialKernel[x_,y_,\[Gamma]_,d_Integer,c_]=(\[Gamma] x.y+c)^d;
(*sigmoid*)
SigmoidKernel[x_,y_,\[Gamma]_,c_]=Tanh[\[Gamma] x.y+c];
(*Vovk's real infinity polynomial universal kernel*)
VovkPolynomialKernel[x_,y_,d_Integer]=(1-x z)^-d;
(*0<q <1 and all compactX \[Subset] ([0, 2 \[Pi] ))^ n*)
FourierKernel[x_,y_,q_]=(1-q^2)/(2 (1-2 q Cos[x-y]+q^2));
(*Wavelet kernel with a \[Element]\[DoubleStruckCapitalR]^ 1 and all compactX\[Subset]\[DoubleStruckCapitalR]^ n*)
WaveletKernel[x_,y_,a_]=Cos[1.75 (x-y)/a] E^(-((x-y)^2/(2 a^2)));

KernelMatrix[K_,X_]:=Outer[K,X,X,1];
KernelMatrix[K_,X_,X1_]:=Outer[K,X,X1,1];



(*different SVM*)
(*the classic SVM, linear separable,C=\[Infinity],
can it be useful?*)
oSVM[X_,Y_,K_,opt___]:=Module[
{aa,ind,ind1,v,w,b,svm},
(*get alpha, linear separable so cs[[2]] is \[Infinity]*)
(*constrains: usually cs[[1]]==0*)
aa=QPSolve[X,Y,K,{\[Infinity],\[Infinity]},opt];
(*Print[aa];*)
ind1=SVindex[aa,Y];
ind=Flatten[Join[ind1]];
(*Print["get v "];*)
v=aa[[ind]]*Y[[ind]];
(*only id k can use w directly*)
If[K==IdentityKernel,
w=v.X[[ind]]];

(*b is calc from all data points, slow..*)
b= -1/2  (Mean[v. KernelMatrix[K,X[[ind]],X[[ind1[[1]]]]]]+ Mean[v.KernelMatrix[K,X[[ind]],X[[ind1[[2]]]]]]);

(*Print[b];*)

svm["offset"]=b;
svm["alpha"]=aa;
svm["svi"]=v;
svm["svX"]=X[[ind]];
svm["svY"]=Y[[ind]];
svm
];

cSVM[X_,Y_,K_,opt___]:=Module[
{aa,ind,ind1,v,w,b,svm,c},
(*C-SVM: cs[[2]] is c*)
c=(cost/.{opt}/.Options[SVMTrain]);
If[c[[1]]<=0\[Or]c[[2]]<=0,Print["option cost error"];Return[$Failed];];
(*constrains: usually cs[[1]]==0*)
aa=QPSolve[X,Y,K,c,opt];
If[aa==$Failed,Return[$Failed]];
(*Print[aa];*)
ind1=SVindex[aa,Y];
ind=Flatten[Join[ind1]];
(*Print["get v "];*)
v=aa[[ind]]*Y[[ind]];

(*b is calc from all data points, slow..*)
b= -1/2  (Mean[v. KernelMatrix[K,X[[ind]],X[[ind1[[1]]]]]]+ Mean[v.KernelMatrix[K,X[[ind]],X[[ind1[[2]]]]]]);

svm["offset"]=b;
svm["alpha"]=aa;
svm["svi"]=v;
svm["svX"]=X[[ind]];
svm["svY"]=Y[[ind]];
svm
];

nuSVM[X_,Y_,K_,opt___]:=Module[
{aa,ind,ind1,v,w,b,svm,n},
n=nu/.{opt}/.Options[SVMTrain];

aa=QPSolve[X,Y,K,n,opt];
If[aa==$Failed,Return[$Failed]];
(*Print[aa];*)
ind1=SVindex[aa,Y];
ind=Flatten[Join[ind1]];
(*Print["get v "];*)
v=aa[[ind]]*Y[[ind]];

(*b is calc from all data points, slow..*)
b= -1/2  (Mean[v. KernelMatrix[K,X[[ind]],X[[ind1[[1]]]]]]+ Mean[v.KernelMatrix[K,X[[ind]],X[[ind1[[2]]]]]]);

svm["offset"]=b;
svm["alpha"]=aa;
svm["svi"]=v;
svm["svX"]=X[[ind]];
svm["svY"]=Y[[ind]];
svm
];



SVindex[\[Alpha]_List]:=
  Flatten[Position[\[Alpha],a_/;a!=0]];
(*this find +1 and -1 sv respectively*)
SVindex[\[Alpha]_List,   y_List]:={Position[\[Alpha]*y,_?Positive]//Flatten,
    Position[\[Alpha]*y,_?Negative]//Flatten};

(*WeightVector[\[Alpha]_List,X_?MatrixQ,y_List]:=
  Sum[y[[i]]*\[Alpha][[ i]]*X[[i]],{i, Length[X]}];*)

(*Bias[\[Alpha]_,X_,Y_,K_]:=
  Module[{sv,ind,v},
 ind=SVindex[\[Alpha],Y];
      -1/2 v. (KernelMatrix[K,X,X[[ind1[[1]]]]]+ KernelMatrix[K,X,X[[ind1[[2]]]]])
];*)


(*this is a wrapper function,
K is kernel, Cs is constrain(corresponding to cost [Cp,Cn]; Cs always >=0, so this is omitted), method is the option for QPsolve algorithm, return $Failed when error occurred*)
QPSolve[X_,Y_,K_,Cs_,opt___]:=Module[
{g,KM,\[Tau]},
(*Print["enter QPSolve"];*)
\[Tau]=t/.{opt}/.Options[SVMTrain];
(*If[0<=\[Tau]<=1,,Return[];];
*)
KM=KernelMatrix[K,X];

Switch[qpType/.{opt}/.Options[SVMTrain],
1,(*classic method*)
QPsolves[KM*Outer[Times,Y,Y],Y,\[Tau],Cs],
2,(*SMO in dlib*)
QPsolvedlib[KM*Outer[Times,Y,Y],Y,\[Tau],Cs],
3,(*special for nu-SVM*)
QPsolvedlib2[KM*Outer[Times,Y,Y],Y,\[Tau],Cs],
_,Print["qpType error."];$Failed
]
];
(*Quadratic problem solution here. Pure Qp
Cs: only Cp used here,
 input Q must be a positive semidefinite matrix, return \[Alpha] *)
(*classical method*)
QPsolves[Q_,Y_,\[Tau]_,cs_]:=Module[
{m,i,j,W,g,vars,sol,p},

m=Length[Y];
(*precision*)
p=Ceiling[-Log[10,\[Tau]]];
If[p<=0,Print["\[Tau] error in QPsolves"];Return[]];

vars=Table[Subscript[\[Alpha], i],{i,1,m}];
(*W=\!\(
\*UnderoverscriptBox[\(\[Sum]\), \(i = 1\), \(m\)]
\*SubscriptBox[\(\[Alpha]\), \(i\)]\)-1/2 \!\(
\*UnderoverscriptBox[\(\[Sum]\), \(i = 1\), \(m\)]\(
\*UnderoverscriptBox[\(\[Sum]\), \(j = 1\), \(m\)]\((Y[\([\)\(i\)\(]\)]\ Y[\([\)\(j\)\(]\)]\ 
\*SubscriptBox[\(\[Alpha]\), \(i\)] 
\*SubscriptBox[\(\[Alpha]\), \(j\)]\ KM[\([\)\(i, j\)\(]\)])\)\)\);*)
W=\!\(
\*UnderoverscriptBox[\(\[Sum]\), \(i = 1\), \(m\)]
\*SubscriptBox[\(\[Alpha]\), \(i\)]\) -1/2 vars.Q.vars;
g=Apply[And,Join[Table[cs[[1]]>= Subscript[\[Alpha], i]>=0,{i,1,m}],{\!\(
\*UnderoverscriptBox[\(\[Sum]\), \(i = 1\), \(m\)]\(Y[\([\)\(i\)\(]\)]\ 
\*SubscriptBox[\(\[Alpha]\), \(i\)]\)\)==0}]];

sol=NMaximize[{W,g},vars,Method->"DifferentialEvolution",AccuracyGoal->p][[2]];

(*m=Max[vars/.sol];*)
(*choose only big enough elements*)
Chop[vars/.sol,\[Tau]]
];
(*this func link to the Dlib lib, since it is fast,\[Tau] can be 0.0001 or smaller *)
(*link;*)
 QPsolvedlib[Q_,Y_,\[Tau]_,cs_]:=Module[
{sol,c,link},
If[!ValueQ[QPsolveC::usage],
link=Install["out`"]];
(*cs's lowb is not used, def 0 here, *)
(*will 2^32 big enough??*)
c={If[cs[[1]]==\[Infinity],N[10^10],N[cs[[1]]]],If[cs[[2]]==\[Infinity],N[10^10],N[cs[[2]]]]};

(*Print[c];*)
(*QPsolved first appears here, so it has context of `PP`*)
sol=QPsolveC[\[Tau],c[[1]],c[[2]],N[Y],Q];
(*Uninstall[link];*)

Chop[sol,\[Tau]]
];
(*this is for nu-SVM,cs has the nu*)
 QPsolvedlib2[Q_,Y_,\[Tau]_,nu_]:=Module[
{sol,c,link},
If[!ValueQ[QPsolveN::usage],
link=Install["out`"]];

(*Print[c];*)
sol=QPsolveN[\[Tau],nu,N[Y],Q];
(*Uninstall[link];*)
(*Print[sol,Options[SVMTrain]];*)
Chop[sol,\[Tau]]
];


End[];
EndPackage[];



